#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include<direct.h>
#include<sys\\stat.h>
#include<sys\\types.h>

#define mSol 1
#define mTerra 3.003167261575587e-6
#define mJupiter 9.543512141169373e-4

#define PI 3.141592653589793
#define G 2.959572553206070e-004
//#define G 1
#define perielioTerra 0.983235919814466
#define vPerielioTerra 6.389661159796174
#define excentricidadeTerra 0.0167

void inicializacao(void);
void artigo2(void);
void criarArquivos(void);
void F3(double *dx, double *z, double *distancia, double *m, int dimensao);
void calcularDistancia(double *z,double *distancia, int dimensao, int corpos);
void jacobiana3C(double* j, double M, double m1,double m2,double r1x,double r1y,double r2x,double r2y,int tamanho);
void produtoEscalarVetor(double* v,double a,int tamanho, double *resultado);
void produtoEscalar(double* v,double* v2,int tamanho, double *resultado);
void somaVetores(double* v,double* v2,int tamanho, double *resultado);
void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado);
void escreverArquivos(void);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);
void finalizar(void);

double *K1,*K2,*K3,*K4;
double  *m,M,tFinal,dt,*y,*Y,*norma,*temp,*temp2;
double  *dx, *distancia,tempo,*J,*z,*z1,*z2,tempo;
int dimensao, n, T, corpos,escrever,tamanho,renormalizacao,orbitas=1,stages=6;
char *exemplo,*aux,*SO;
clock_t inicio,fim,total;
FILE *f[6];

int main(int argc, char *argv[])
{
	inicio=clock();

	escrever=atoi(argv[1]);
	tFinal=atof(argv[2]);
	n=atof(argv[3]);

	inicializacao();

	calcularDistancia(z,distancia,dimensao,corpos);
	escreverArquivos();

	//RK4
	int i,j;
	for(i=1;i<T;i++)
	{
		//temp2=soma parcial

		//RK4
		//stage 1
		F3(temp, z, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K1);

		//stage 2
		produtoEscalarVetor(K1,0.5,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K2);

		//stage 3
		produtoEscalarVetor(K2,0.5,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K3);

		//stage 4
		produtoEscalarVetor(K3,1.0,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F3(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K4);

		produtoEscalarVetor(K1,1.0/6,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K2,1.0/3,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K3,1.0/3,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K4,1.0/6,tamanho,temp);
		somaVetores(z,temp,tamanho,z);

		norma[1]=calculaNorma(z+tamanho/2,tamanho/2);

		y[1]=((double)(i-1)/(double)(i))*y[0]+2*log(norma[1]/norma[0]);
		Y[1]=((double)(i-1)*Y[0]+y[1])/(double)i;

		if(i%renormalizacao==0)
		{
			for(j=tamanho/2;j<tamanho;j++)
			{
				z[j]=z[j]/norma[1];
			}
			norma[1]=calculaNorma(z+tamanho/2,tamanho/2);
		}

		if(i%escrever==0)
		{
			escreverArquivos();
		}

		y[0]=y[1];
		Y[0]=Y[1];
		norma[0]=norma[1];
	}

	printf("\n%lf\n",Y[1]);

	int minutos=0;
	int segundos=0;
	int horas=0;
	tempo=0;
	int tempo2=0;

	fim=clock();
	total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);

	sprintf(aux,"%0.10lf ",tempo);
	fputs(aux,f[2*corpos]);
	fputs("\n",f[2*corpos]);

	finalizar();

	return 0;
}

void inicializacao(void)
{
	dimensao=2;
	corpos=2;
	renormalizacao=100;
	T=(int)pow(2.0,n)+1;
	dt=tFinal/(T-1);
	tamanho=corpos*dimensao*2*2;

	int i;

	m =(double *)calloc(corpos, sizeof(double));
	z =(double *)calloc(tamanho, sizeof(double));
	temp =(double *)calloc(tamanho, sizeof(double));
	temp2 =(double *)calloc(tamanho, sizeof(double));
	dx =(double *)calloc(tamanho, sizeof(double));
	distancia =(double *)calloc(corpos*corpos, sizeof(double));
	Y=(double*)calloc(2,sizeof(double));
	y=(double*)calloc(2,sizeof(double));
	norma=(double*)calloc(2,sizeof(double));
	J=(double *)calloc(pow(tamanho/2,2), sizeof(double));
	K1=(double*)calloc(tamanho,sizeof(double));
	K2=(double*)calloc(tamanho,sizeof(double));
	K3=(double*)calloc(tamanho,sizeof(double));
	K4=(double*)calloc(tamanho,sizeof(double));
	aux =(char *)calloc(1000, sizeof(char));

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=tamanho/2;i<tamanho;i++)
	{
		z[i]=rand()/(double)RAND_MAX;
		z[i]=z[i]-0.5;
		//z[i]=1.0;
	}

	// norma delta inicial nao normalizado
	norma[0]=calculaNorma(z+tamanho/2,tamanho/2);
	//normalizacao delta inicial
	for(i=tamanho/2;i<tamanho;i++)
	{
		z[i]=z[i]/norma[0];
	}
	imprimir(z+tamanho/2,tamanho/2);
	//norma delta inicial normalizado
	norma[0]=calculaNorma(z+tamanho/2,tamanho/2);

	//somaVetores(z,z+2*dimensao*corpos,2*dimensao*corpos,z);

	artigo2();
	criarArquivos();
}

void artigo2()
{
	exemplo="C:\\dados\\";
	mkdir(exemplo);
	exemplo="C:\\dados\\artigo2EP";
	mkdir(exemplo);

	//corpo0
	M=0.3*mSol;

	//corpo1

	m[0]=mJupiter;
	double periodo=7.5;
	z[0*dimensao+0]=G*M*pow(periodo,2)/(4*pow(PI,2));
	z[0*dimensao+0]=pow(z[0*dimensao+0],((double)1.0/3));
	z[2*dimensao+1]=m[0]*sqrt(G*M/z[0*dimensao+0]);

	//corpo2
	m[1]=mTerra;

	//posicao inicial corpo 2:
	//a2=semi-eixo maior, b=semi-eixo menor, c= distância focal,
	double a=0.05;
	double e=0.7;
	double c=e*a;
	double b=sqrt(pow(a,2)-pow(c,2));
	z[1*dimensao+0]=a*(1-e);

	double P=sqrt(4*pow(PI,2)*pow(a,3)/(G*M));
	double N=2*PI/P;
	z[3*dimensao+1]=-1*m[1]*a*b*N/z[1*dimensao+0];
}

void criarArquivos(void)
{
	char *raiz,*diretorio,*aux,*barra;
	int i;

	raiz=(char*)calloc(100,sizeof(char));
	aux=(char*)calloc(100,sizeof(char));
	diretorio=(char*)calloc(100,sizeof(char));

	barra="\\";
	strcpy(raiz,exemplo);
	strcat(raiz,barra);
	mkdir(raiz);
	sprintf(aux,"%0.0lf",tFinal);
	strcat(raiz,aux);
	strcat(raiz,barra);
	mkdir(raiz);
	sprintf(aux,"%d",n);
	strcat(raiz,aux);
	strcat(raiz,barra);
	mkdir(raiz);

	if(orbitas)
	{
		for(i=0;i<corpos;i++)
		{
			strcpy(diretorio,raiz);
			strcat(diretorio,"r");
			sprintf(aux,"%d",(i+1));
			strcat(diretorio,aux);
			strcat(diretorio,".txt");
			f[i] = fopen(diretorio,"w");
			if(f[i] == NULL)
			{
				printf("Erro ao abrir arquivo!!!\n");
				exit(1);
			}
			strcpy(diretorio,raiz);
			strcat(diretorio,"p");
			sprintf(aux,"%d",(i+1));
			strcat(diretorio,aux);
			strcat(diretorio,".txt");
			f[i+corpos] = fopen(diretorio,"w");
			if(f[i+corpos] == NULL)
			{
				printf("Erro ao abrir arquivo!!!\n");
				exit(1);
			}
		}
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"tempo");
	strcat(diretorio,".txt");
	f[2*corpos] = fopen(diretorio,"w");
	if(f[2*corpos] == NULL)
	{
		printf("Erro ao abrir arquivo5!!!\n");
		exit(1);
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"megno");
	strcat(diretorio,".txt");
	f[2*corpos+1] = fopen(diretorio,"w");
	if(f[2*corpos+1] == NULL)
	{
		printf("Erro ao abrir arquivo megno!!!\n");
		exit(1);
	}

	free(aux);
	free(raiz);
	free(diretorio);

}

void F3(double *dx, double *z, double *distancia, double *m, int dimensao)
{

	dx[0]=z[4]/m[0];
	dx[1]=z[5]/m[0];
	dx[2]=z[6]/m[1];
	dx[3]=z[7]/m[1];

	calcularDistancia(z,distancia,dimensao,corpos);

	dx[4]=      (G*m[0]*M*(-z[0]))/(double)(pow(sqrt(z[0]*z[0]+z[1]*z[1]),3));
	dx[4]=dx[4]+(G*m[0]*m[1]*(z[2]-z[0]))/(double)(pow(distancia[0*(corpos)+1],3));
	dx[5]=      (G*m[0]*M*(-z[1]))/(double)(pow(sqrt(z[0]*z[0]+z[1]*z[1]),3));
	dx[5]=dx[5]+(G*m[0]*m[1]*(z[3]-z[1]))/(double)(pow(distancia[0*(corpos)+1],3));

	dx[6]=       (G*m[1]*M*(-z[2]))/(double)(pow(sqrt(z[2]*z[2]+z[3]*z[3]),3));
	dx[6]=dx[6]+(G*m[1]*m[0]*(z[0]-z[2]))/(double)(pow(distancia[1*(corpos)+0],3));
	dx[7]=       (G*m[1]*M*(-z[3]))/(double)(pow(sqrt(z[2]*z[2]+z[3]*z[3]),3));
	dx[7]=dx[7]+(G*m[1]*m[0]*(z[1]-z[3]))/(double)(pow(distancia[1*(corpos)+0],3));

	jacobiana3C(J,M,m[0],m[1],z[0], z[1], z[2], z[3],2*dimensao*corpos);

	produtoMatrizVetor(J,z+8,2*corpos*dimensao,dx+8);
}

void calcularDistancia(double *z,double *distancia, int dimensao, int corpos)
{
	int i,j,k;
	double temp;

	for(i=0;i<corpos;i++)
	{
		for(j=0;j<corpos;j++)
		{
			distancia[i*(corpos)+j]=0;
			for(k=0;k<dimensao;k++)
			{
				temp=(z[i*dimensao+k]-z[j*dimensao+k]);
				temp=temp*temp;
				distancia[i*(corpos)+j]=distancia[i*(corpos)+j]+temp;
			}
			distancia[i*(corpos)+j]=(double)sqrt((double)distancia[i*(corpos)+j]);
		}
	}
}

void jacobiana3C(double* j, double M, double m1,double m2,double r1x,double r1y,double r2x,double r2y,int tamanho)
{
	double t2 = 1.0/m1;
	double t3 = 1.0/m2;
	double t4 = r1x-r2x;
	double t5 = r1y-r2y;
	double t6 = r1x*r1x;
	double t7 = r1y*r1y;
	double t8 = t6+t7;
	double t13 = r1x*2.0;
	double t14 = r2x*2.0;
	double t9 = t13-t14;
	double t10 = t4*t4;
	double t11 = t5*t5;
	double t12 = t10+t11;
	double t15 = 1.0/pow(t12,(5.0/2.0));
	double t16 = 1.0/pow(t8,(5.0/2.0));
	double t17 = 1.0/pow(t12,(3.0/2.0));
	double t18 = t9*t9;
	double t19 = G*m1*m2*t15*t18*(3.0/4.0);
	double t20 = r1y*2.0;
	double t21 = r2y*2.0;
	double t22 = t20-t21;
	double t23 = G*m1*m2*t15*t22*t9*(3.0/4.0);
	double t24 = G*M*m1*r1x*r1y*t16*3.0;
	double t25 = t23+t24;
	double t26 = G*m1*m2*t17;
	double t27 = 1.0/pow(t8,(3.0/2.0));
	double t28 = t22*t22;
	double t29 = G*m1*m2*t15*t28*(3.0/4.0);
	double t30 = -t19+t26;
	double t31 = r2x*r2x;
	double t32 = r2y*r2y;
	double t33 = t31+t32;
	double t34 = 1.0/pow(t33,(5.0/2.0));
	double t35 = t26-t29;
	double t36 = G*M*m2*r2x*r2y*t34*3.0;
	double t37 = t23+t36;
	double t38 = 1.0/pow(t33,(3.0/2.0));

	int colunas=2*2*2;

	j[0*colunas+4]=1/m1;
	j[1*colunas+5]=1/m1;
	j[2*colunas+6]=1/m2;
	j[3*colunas+7]=1/m2;

	j[4*colunas+0]=t19-G*M*m1*t27-G*m1*m2*t17+G*M*m1*t16*t6*3.0;
	j[4*colunas+1]=t25;
	j[4*colunas+2]=t30;
	j[4*colunas+3]=-t23;

	j[5*colunas+0]=t25;
	j[5*colunas+1]=-t26+t29-G*M*m1*t27+G*M*m1*t16*t7*3.0;
	j[5*colunas+2]=-t23;
	j[5*colunas+3]=t35;

	j[6*colunas+0]=t30;
	j[6*colunas+1]=-t23;
	j[6*colunas+2]=t19-t26-G*M*m2*t38+G*M*m2*t31*t34*3.0;
	j[6*colunas+3]=t37;

	j[7*colunas+0]=-t23;
	j[7*colunas+1]=t35;
	j[7*colunas+2]=t37;
	j[7*colunas+3]=-t26+t29-G*M*m2*t38+G*M*m2*t32*t34*3.0;

}

void produtoEscalar(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	*resultado=0;
	for(i=0;i<tamanho;i++)
	{
		*resultado+=v[i]*v2[i];
	}
}

void somaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]+v2[i];
	}
}

void produtoEscalarVetor(double* v,double a,int tamanho,double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=a*v[i];
	}
}

void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado)
{
	int i,j;
	double temp;

	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		resultado[i]=temp;
	}
}

void escreverArquivos(void)
{
	int j,k;
	char aux[100];

	if(orbitas)
	{
		for(k=0;k<corpos;k++)
		{
			for(j=0;j<dimensao;j++)
			{
				//arquivo txt
				//posicao
				sprintf(aux,"%0.20lf ",z[k*dimensao+j]);
				fputs(aux,f[k]);
				//momento
				sprintf(aux,"%0.20lf ",z[(k+corpos)*dimensao+j]);
				fputs(aux,f[k+corpos]);
			}
			fputs("\n",f[k]);
			fputs("\n",f[k+corpos]);
		}
	}
	//megno
	sprintf(aux,"%0.15lf ",Y[1]);
	fputs(aux,f[2*corpos+1]);
	fputs("\n",f[2*corpos+1]);

}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}

void finalizar(void)
{
	if(orbitas)
	{
		fclose(f[0]);
		fclose(f[1]);
		fclose(f[2]);
		fclose(f[3]);
	}
	fclose(f[4]);
	fclose(f[5]);

	free(distancia);
	free(m);
	free(K1);
	free(K2);
	free(K3);
	free(K4);
	free(J);
	free(temp);
	free(norma);
	free(Y);
	free(y);
	free(z);
}
